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Abstract 

When an unmagnetized plasma comes in contact with a material surface, the difference in mobil¬ 
ity between the electrons and the ions creates a non-neutral layer known as the Debye sheath (DS). 
However, in magnetic fusion devices, the open magnetic field lines intersect the structural elements 
of the device with near grazing incidence angles. The magnetic field tends to align the particle flow 
along its own field lines, thus counteracting the mechanism that leads to the formation of the DS. 
Recent work using a fluid model [P. Stangeby, Nucl. Fusion 52 , 083012 (2012)] showed that the 
DS disappears when the incidence angle is smaller than a critical value (around 5° for ITER-like 
parameters). Here, we study this transition by means of numerical simulations of a kinetic model 
both in the collisionless and weakly collisional regimes. We show that the main features observed 
in the fluid model are preserved: for grazing incidence, the space charge density near the wall is 
reduced or suppressed, the ion flow velocity is subsonic, and the electric field and plasma density 
profiles are spread out over several ion Larmor radii instead of a few Debye lengths as in the unmag¬ 
netized case. As there is no singularity at the DS entrance in the kinetic model, this phenomenon 
depends smoothly on the magnetic field incidence angle and no particular critical angle arises. The 
simulation results and the predictions of the fluid model are in good agreement, although some 
discrepancies subsist, mainly due to the assumptions of isothermal closure and diagonality of the 
pressure tensor in the fluid model. 
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I. INTRODUCTION 


In magnetic fusion devices such as tokamaks, the conhning magnetic held is designed 
so that the held lines that intersect some machines components do so with near grazing 
incidence in order to maintain power deposition within sustainable limits. Due to the large 
diherence in inertia between the ions and the electrons, the latter tend to be lost to the 
absorbing wall faster than the former, leading to the formation of a thin (a few Debye 
lengths wide) positively-charged transition layer in front of the wall, the so-called Debye 
sheath (DS) (see [1] for a large-scope review on the topic). The resulting large electric held 
in the DS repels the electrons and accelerates the ions, leading to a sustainable steady-state 
with zero net current at the wall. 

In the presence of a magnetic held whose direction is not normal to the wall, the structure 
of the transition is more intricate. The magnetic held maintains the ions how aligned with 
its own direction, while the electric held tends to accelerate them normally to the wall, 
leading to a competition between these two ehects. In the case of nearly grazing incidence, 
the particle motion along the normal to the wall is essentially cyclotronic, resulting in a 
strongly reduced net how in that direction. The efficiency of the conhnement decreases 
when one approaches the wall, as more and more Larmor orbits intersect the wall. As the 
electrons are more strongly conhned than the ions, there exists a new transition layer, the 
so-called Chodura sheath (CS) or magnetic pre-sheath [2], where the imbalance between 
the ionic and electronic hows is sufficiently compensated by the diherence in conhnement to 
maintain quasi-neutrality. This transition layer, between a fully magnetized plasma how and 
the wall is typically a few ion Larmor radii thick. Since generally p* S> Ad the plasma-wall 
transition is globally smoother than in the purely electrostatic case, with smaller spatial 
gradients for the electric held and plasma density near the wall. 

In the most general case, the DS and the CS coexist: the imbalance between the ionic and 
electronic parallel how still requires the formation of a positively charged DS in order to 
ensure ambipolarity at the wall. The boundary between the CS and the DS is characterized 
by the breakdown of quasi-neutrality and the onset of a supersonic ion how velocity at 
the entrance of the DS. For unmagnetized plasmas, this reduces to the well-known Bohm 
criterion PS]. A similar criterion was derived by Chodura [2] in the magnetized case, which 
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requires the parallel ion flow velocity at the entrance of the CS to be supersonic. 

In the landmark study by Chodura [2], the main features of the CS-DS transition were 
described using both a fluid model and numerical results from particle-in-cell (PIC) simula¬ 
tions. Further studies of the plasma-wall transition, focussing on its stability, were performed 
with PIC simulations [SIE]. The fluid model was later extended with friction terms to en¬ 
compass both the magnetic and collisional presheath [7] (and more recently [8]). This model 
was recently used to show some partial agreement with experimental data ^ in a different 
regime {\coii ^ Pi ^ ^d) with respect to the one considered here {\coii ^ Pi ^ Xd)- 

In a recent work [lOj, Stangeby also used a fluid model to examine the CS-DS transition 
for low values (a few degrees) of the incidence angle of the magnetic held, ie, in the range 
relevant to the plasma-divertor interaction in fusion devices. Importantly, this study showed 
the existence of a critical incidence angle under which the plasma-wall transition occurs 
without the need for the formation of the DS. As a result, the electric held and the plasma 
density gradients are not restricted to the (very thin) DS, but extend much further (a few ion 
Larmor radii) into the CS. This effect is significant enough to have a non-negligible impact 
on prompt redeposition of sputtered neutrals in a tokamak scrape-off-layer (SOL) [TTl 112] . 


This potentially important application warrants a more detailed analysis of this phe¬ 
nomenon, going beyond the simple fluid approach that was used in Ref. [10]. The main 
objective of the present paper is to examine the robustness of Stangeby’s results by means 
of numerical simulations of a kinetic model [13]. Various effects that can have an impact 
on the transition will be analyzed in details, such as the magnitude and incidence of the 
magnetic held, the effect of collisions, and isotopic effects. Generally speaking, Stangeby’s 
results are confirmed: the DS disappears for small angles of incidence (1° — 5°), although 
the transition is not as clear-cut as in the fluid model. Note that we will not consider here 
the extreme case a < rrielmi 1° (for deuterium), for which the ions reach the wall faster 
than the electrons, and consequently the sheath structure changes considerably [IT] . 


The present paper is organized as follows: In Sect. |TT[ we summarize the results obtained 


by Stangeby using a fluid model. In Sect. [ITT| we describe the kinetic model and the 
numerical method and parameters. In Sect. |IV[ we examine the CS-DS transition using a 
collisionless model, with parameters and boundary conditions chosen to match as closely as 
possible those of Ref. [lO]. In Sect. |V| we directly compare the spatial profiles obtained 
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Figure 1: Geometry of the model 


from the fluid model and kinetic simulations. In Sect. |VI[ we introduce a collision operator 
in our kinetic model, and use it to check the robustness of the observations made in the 


collisionless regime. In section VII, we summarize the main conclusions of this study and 
mention some of the key issues that remain to be addressed. 


II. STANGEBY’S RESULT FROM FLUID THEORY 

Stangeby [10] considered a plasma composed of electrons of charge —e and a single ion 
species of charge g* = ZjC. The plasma is bounded by a fully absorbing wall on one side, 
while thermal equilibrium is assumed far from the wall (see Fig. [^. Noting x the direction 
corresponding to the normal to the wall, the system is assumed invariant by translation in 
the (g, z) plane parallel to the wall. The plasma is magnetized by an external magnetic field 
Bq, constant in space and time, whose direction is normal to and makes an angle a with 
the wall, i.e Bq = i?o(sinae 3 , + cosae^). The self-consistent magnetic field generated by 
plasma currents is neglected. 

The main result of Ref. [10] is the existence of a critical angle Oc for which there is strictly 
no Debye sheath, or more precisely the average flow along the normal to the wall never 
becomes sonic. We will first reestablish this result with slightly more relaxed assumptions 
in order to treat both sonic and supersonic regimes, and then examine the actual simulation 
results. Though the model used in is a fluid one, the result is actually quite generic. 
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From the ion flux conservation dxjxi = 0 we have : 


{Vx)i{x) = 


■w 

J XI 


Zij, 


w 

xi 


ni{x) ne{x)+p/e 


( 1 ) 


where {vx)i is the mean ion velocity, jxi is the mean ion current, is the ion (electron) 
number density, and p = e{Zini — Ue) is the charge density. The superscript “W” refers to 
the wall and (•) stands for the averaging operator over velocity space. Using the ambipolarity 
condition at the wall we have 

w 

T? 

{vx)i{x) = [(n||)f sin a + (v^ ■ cos a] +p{x)/e ' 


Now we make two assumptions. The first is on the ratio — / / u , which we will take to 

be less than unity. Such condition is fulfilled in the case of a quasi-neutral region (p cs 0, 
as in the CS) or a positively charged region (p > 0, as in the standard DS), subject to the 
condition of a decrease of the electron density when one approaches the wall {dxUe < 0). This 
is clearly the case for Boltzmann electrons and a negatively charged wall, as was assumed 
in Ref. [10]. Whatever the exact assumptions, as long as „ {x)+p{x)/e ^ 1 we obtain a bound 
on the ion flow velocity 


|(w)i(a:)| < |(t'||)r sin a (v_l ■ e^,)^ cos a 


( 3 ) 


The second assumption is that the electrons are perfectly magnetized up to the wall, ie, 
(v_L ■ = 0. This becomes obviously false for distances smaller than the electron Larmor 

radius pe from the wall, but can be considered a reasonable approximation as long as the 
electron flow variation is mild. We then have 


|(nr.)i(a:)| < sina|(n||)f 


( 4 ) 


For sufficiently small a, the bound of Eq. (|^ may prevent the ion mean velocity (vx)i 
from becoming supersonic, in which case no DS is required to guarantee ambipolarity. This 
happens when a is equal or smaller than the critical value ac defined as 


sincTc 



( 5 ) 


In the case of a half-Maxwellian electron parallel velocity distribution at the wall, one has 
('^11 = y/Teo/ (27ime) and the result of Ref. [10] is readily obtained. The underlying 
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physical phenomenon is essentially the limitation of the electron current at the wall by 
the magnetic held, which entails a limitation of the ion current. For sufficiently small 
a, an ambipolar how along x can be maintained at the wall without requiring strong ion 
acceleration, so that there is no need for a DS. 

A few points of importance should be noted: 

1. While the bound on the CS ion how velocity in Eq. (|^ is quite generic, the notion of 
a well-dehned critical angle stems from two assumptions: a Bohm criterion on the ion 
velocity for the existence of the sheath (ie, |(n 3 ;)i| > Cs at the sheath entrance) and 
perfect magnetization of the electrons. In a kinetic model such as the one considered 
later on in this paper, the relationship between the mean ion how velocity and the 
sheath stability is not as direct as the simple Bohm criterion. 

2. A second point is the fact that the bound of Eq. Q and the critical angle do not 
depend explicitly on the how at the CS entrance, and are thus valid in the CS in 
both the sonic and supersonic regimes. This is in contrast with the result presented 
in Appendix A of Ref. [lO] which relies on the erroneous use in a supersonic case of 
the potential drop in the CS that had been established for a sonic case (Eq. (33) in 
[Tn] . used in conjunction with Eq. (A3) of the same paper). 

3. As was noted in [10], in a model accounting for the hnite electron Larmor radius, the 
angular dependency of the electron current would be more complex than the simple 
sin a behaviour considered here. 

III. KINETIC MODEL AND NUMERICAL PARAMETERS 

In the kinetic model considered here, the dynamics of the ions is described by the evolu¬ 
tion of the phase-space distribution function fi(t, x,Vx,Vy,Vz) obeying the collisional Vlasov 
equation 




where Ud = ZieBo/rrii is the ion cyclotron frequency. In all results presented hereafter 
the collision operator, whenever present, is a Bathnagar-Krook-Gross (BGK) linear re¬ 
laxation operator |T5], which drives the distribution function to an isotropic Maxwellian 
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distribution, ie, Ci{fi) = —r'i(/ — f^) where Ui is the ion relaxation rate and = 
/ \3/2 r 2-1 

,n I 1 exp — . At the wall (x = 0), an absorbing boundary condition is as- 


^iO 


\ 27rTio 




sumed in x for the incoming part of the distribution function, ie, fi{t,0,Vx,Vy,Vz) = 0 for 
Vx > 0. On the plasma side {x = L), the incoming particle distribution is prescribed by 
fi{t, L,Vx,Vy,Vz) = f-"'{vx,Vy,Vz) for Vx < 0. In the collisional simulations /*” is simply a 
Maxwellian with bulk plasma parameters (the same that is used for the BGK operator). In 
the collisionless simulations, it is a held-ahgned drifting distribution with parallel velocity 


that satisfies the Chodura criterion at the CS entrance (see Sect. IVA). 


The electrostatic held E = —dx4> is computed from the electrostatic potential by solving 
the Poisson equation 

dlx +-iZirii-rie) =0 ( 7 ) 

^0 

with a Dirichlet boundary condition (p = 0 at x = L and a Von Neumann condition Ex = cr/eo 
at the wall. The wall charge surface a is computed by integrating in time the outgoing net 
electric current: 

t 

cr = -e y ^ Zsjxs{t\x = Q)dt\ 

Q s=i^e 

with jxs = / Vxfsd^v. The full electron kinetic dynamics is not resolved, but instead a 
Boltzmann law is assumed for the electron density Ug = u^e/exp [e(0 — 0j.e/)/77o)]- The 
reference quantities are dehned at x = L by (pref = 0 and riref = ni{L). The outgoing 
electron hux at the wall is computed by assuming a half-Maxwellian distribution and is 
given by 

jZ = - sina^exp [e(0 - (pref)/Teo] . (8) 

The latter relation does not take into account hnite electron Larmor radius effects, as it is 
assumed that jxe = sinajue. 


All numerical simulations were performed using the Eulerian code described in Ref. [T3] . 
The numerical scheme is based on a split-operator technique for the time-stepping algorithm, 
with interpolations performed with a positive hux conservative (PEG) scheme [I6]. In all 
cases, starting from a uniform Maxwellian plasma, the system is left to relax towards a 
stationary state. A hrst set of simulations were run in a collisionless regime (z/j = 0) over 
a spatial domain limited to the GS-I-DS region, covering a few ion Larmor radii. A second 
set of simulations were run in a collisional regime where the full transition from an isotropic 
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Maxwellian plasma to the wall is considered, including the collisional presheath. In both 
cases, parametric scans with a G {2°, 3°, 4°, 5°, 10°, 15°, 30°, 45°, 60°, 90°} were performed. 


IV. COLLISIONLESS PLASMA-WALL TRANSITION 

The parameters of the hrst set of simulations were set in order to match as closely as 
possible those of the fluid model used in Ref. [lO]. The simulation box length is between 
L K, 120Ad and L 800Ad, depending on the strength of the magnetic held (in Stangeby’s 
quasi-neutral model, since the Debye length vanishes, the CS entrance is rejected at inhnity). 
Parametric scans in the incidence angle a were performed for hydrogen {rrii = mu) and 
deuterium {rrii = ^mu)- The magnetic held intensity is such that = O.ObWpj, ie pi = 20Ad- 
For all simulations, we assumed equal temperatures Tjo = T^q. For brevity, the local value 
of any quantity X expressed at the wall {x = 0) and at the magnetic presheath entrance 
{x = L) will be tagged respectively as and X^^^. 


A. Boundary conditions 


At the plasma boundary, ie, the CS entrance, the incoming ion hux should be supersonic 
(Chodura criterion) and aligned with the magnetic held direction. To this end, we prescribe 
the following distribution function at x = L\ 


/r = M(-^||) ^ exp 


■> (S!)' 


2 ?^ 


(9) 


where H is the Heaviside function, vm = y/Tio/rrii, K = ^yp-2 2 ^F(^^), with F the 
Euler gamma function. The average parallel velocity corresponding to fl'^ is (f||) = 
—nt/ii\/2F(^^)/F(^^). In the results presented in this section the (5 exponent was set 
equal to 2, leading to an average how (n||) = —l.Qvthi) ie, slightly supersonic. Smaller values 
of the parallel velocity may run the risk of destabilizing the transition. The above distribu¬ 
tion is compatible with huid models that assume a sonic or slightly supersonic how at the 
entrance of the CS. In addition, it is also compatible with the velocity distribution obtained 
self-consistently from a kinetic model that incorporates weak collisions, as we shall show in 
Sect. El 
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(a) (b) 

Figure 2: (a) Prescribed ion distribution functions at x = L averaged over {vy, Vz), for a few values 
of the incidence angle a (the distributions are normalized to the peak value of the a = 90° case); 
(b) Total potential drop in the CS and DS as a function of a. 


In Fig. the Vx dependency of the incoming distribution function is shown for a few values 
of a. The case a = 90° corresponds to Vx = f||. One should note that the parallel velocity 
distribution is not a Maxwellian, and that its effective ’’temperature” T*”'' = ~ 0.45Tjo 

is smaller than T^q. In a magnetic-held-aligned basis such as (bje^ x b,e^), the kinetic 
pressure tensor is diagonal but anisotropic. In the {x, y, z) basis, it is not even diagonal 
anymore and the various components of the pressure tensor vary with a. For instance, the 
XX component of the pressure tensor, for = 2, is equal to ~ njoTio(l — 0.55 sin^ a). 

In a collisionless model, the total potential drop from the CS entrance to the wall is 
independent on the angle a. However, for very small angles, numerical errors (due to the 
presence of a small but non-zero electric held near x = L) slightly break this invariance. 
This entails a small variation with a of the total potential drop, as seen in Fig. [^. However, 
this small error does not affect the main physical conclusions that can be drawn from the 
forthcoming numerical simulations. 
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(a) 

Figure 3: Spatial profiles of the charge density (. 
case with deuterium ions. 



(b) 

.) and the electric potential (b), for a collisionless 




(a) (b) 

Figure 4: Spatial prohles of the electrostatic field (a) and the ion density (b), for a collisionless 
case with deuterium ions. 


B. Effect of the angle of incidence 


We will now consider the parametric dependency of the CS-DS transition with the magnetic 
held incidence angle a. Figurej^ shows that the space charge density near the wall decreases 
rapidly with decreasing a. Although the charge density does not strictly vanish (nor changes 
sign), the strong limitation of the space charge density is a clear signature that the DS 
progressively disappears at small incidence angles. In addition, the spatial prohle of the 
electric potential (Fig. [^) evolves from a two-scale prohle at large a. - typical of the CS-DS 
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Figure 5: Spatial profiles of the average ion flow velocity perpendicular to the wall, for a collisionless 
case with deuterium ions. 

transition - to a smooth evolntion at low a. As a conseqnence, althongh the peak of the 
electric held decreases strongly as the DS vanishes (Fig. |^), its extension reaches mnch 
farther into the plasma, several ion Larmor radii from the wall. As discnssed in Ref. [TO] , 
this is of signihcant importance for the estimation of the prompt redeposition of spattered 
imparity ions: indeed, while the overall electric held intensity decreases with a, it will ahect 
spattered nentrals ionized farther from the wall and increase prompt redeposition. 

The ion (and thus plasma) density drop is also spread out and reaches lower values with 
decreasing a (Fig. |^). This depletion of the plasma density near the wall (for regions such 
that X < Pi) entails a lower ionization rate for sputtered neutrals, thus lowering prompt 
redeposition. 

Let us now consider the ion mean how perpendicular to the wall (Fig. [^. Due to both the 
anisotropic nature of the kinetic pressure tensor and the non-uniformity of the ’’tempera¬ 
tures” (see Sect. |V]for a discussion of the huid closure), we refrain here from normalizing the 
how to the usual sound speed Cg = y/(T^o -|- Tf,Q)/nii ~ lAvthi, which is strictly valid only 
in the case of an isothermal closure for the P^x component of the kinetic pressure tensor. 
In our case, the sound speed can be roughly estimated (from f-^) as ranging from l.2vthi 
to lAvthi when a ranges from 90° to 2°, and is very close to lAvthi for the lowest range 
(a < 15°) of angles considered. 

Figure clearly shows that the peak (absolute) value of the ion mean velocity decreases 
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(a)a = 3° 


(b)a = 30° 



(c)a = 90° 

Figure 6: Ion velocity distribution function in Vx for several positions (indicated on top of each 
curve, in units of Xd) and three values of a: (a) a = 3°, (b) a = 30°, and (c) a = 90°. For 
each value of a, all distributions are normalized to their peak value at the entrance of the CS 

(x = 117.3 Ad). 


with decreasing a and is limited to snbsonic values for low angles of incidence, below 
approximately 5°. Together with the disappearance of the space charge in front of the wall 
(Fig. these results confirm Stangeby’s conclusion that no DS forms below a certain 
angle of incidence. The limitation of both the ion density and the average velocity with 
decreasing a are clearly visible when examining directly the Vx velocity profile of the ion 
distribution function (averaged over Vy and v^), as shown in Fig. 

We will now examine more closely the behaviour with a of a few important quantities 
measured at the wall. The x component of the electrostatic held at the wall is shown in 
Fig. 0 as a function of sin a. As expected from the above observations, it is an increasing 
function of a. For the smallest angles a G {2°, 3°, 4°, 5°, 10°} (inset of Fig. [^, the evolution 
is roughly linear in sin a, but the overall behaviour for the full range of angles is less obvious. 

The space charge density at the wall clearly exhibits a linear dependency in sin a (Fig. [^). 
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Figure 7: Electric field magnitude at the wall as a function of the incidence angle, for a collisionless 
case with hydrogen or deuterium ions. The inset is a zoom at small angles. 


This fact allows us to obtain a semi-empirical fit for the ion perpendicular velocity at the 
wall. Indeed, taking Eq. ([^ at the wall with an electron current = — sin a{vthe/ 
we obtain 


\{v.)Y\ = 


Vthe 


sma 


Vthe 


sm a 


( 10 ) 


\/27r 1 + + 

where k is a fitting parameter. To obtain Eq. we have assumed that cx sina (see 

Fig. [^) and that is independent of a. An interesting fact here is that the coefficient 
K can be computed in the normal incidence case (a = 90°), which does not require a full 
1D3V model but only a far simpler IDIV simulation. Once k has been determined, the ion 


perpendicular flow for any incidence angle can be computed using Eq. (10). Some examples 
of this semi-empirical fit are shown in Fig. if. for both hydrogen and deuterium ions. 


C. Effect of the magnetic field amplitude at fixed angle (a = 2°) 

In the simulations considered so far, the scaling pi = 20Xd 3> \d (or cOci/cOpi = 0.05 1) 

was valid. In that regime, decreasing the magnetic field intensity Bq will essentially result 
in a rescaling of the CS, whose thickness increases with growing pi. This is clear from Fig. 

for instance in the case ujd = O.OlWpj (p* = lOOAij), where the velocity profile stretches 
out to several hundred Debye lengths. At the same time, the charge separation near the 
wall decreases with decreasing magnetic field, and almost disappears for uJd = O.OloOpi (Fig. 
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(a) (b) 

Figure 8: Charge density (a) and ion mean velocity (b) at the wall as a function of the incidence 


angle. The coefficient k of Eq. (10) is obtained from the simulation data for a = 90°. 


10). This is because, the CS being larger, the whole potential drop can more easily occur 


within the CS, with hardly any need for a non-neutral DS. Thus, in the large Pi/Xo regime, 
the disappearance of the DS predicted by Stangeby is even more apparent. 

In contrast, increasing Sq, and thus decreasing pi, results in a progressive breaking of the 
above scaling (see Ref. im for a discussion of the scales entering the transition). For the 
case of low incidence angles, the consequences are twofold. On the one hand, we observe 
a stronger limitation of the ion flux perpendicular to the wall, as can be seen from Fig. 
1^ On the other hand, the charge separation near the wall tends to increase with Bq (Fig. 


10). These observations can be explained as follows. With increasing Bq, the CS extension 


becomes of the same order as that of the DS, so that the two sheaths overlap. Since the 
total potential drop remains constant, the overall width of the transition zone becomes too 
narrow to allow a quasi-neutral transition. Consequently, the almost quasi-neutral transition 
previously observed for low magnetic fields at grazing incidence (curve corresponding to 


^ci/^pi = 0.05 in Fig. 10) disappears, and the formation of a sheath is again required to 
ensure a smooth plasma-wall transition. This effect may be interpreted as the appearance 
of a “new” type of non-neutral sheath, whose thickness is of the order of the ion Larmor 
radius, when the scaling pi ^ Xjy is satished. 


14 


















































Figure 9: Ion flow velocity perpendicular to the wall for various amplitudes of the magnetic held 
-Bo and a = 2°. Collisionless simulations with deuterium ions. 



Figure 10: Charge density profile for various amplitudes of the magnetic field Bq and a = 2°. 
Collisionless simulations with deuterium ions. 


D. Non-floating (biased) wall 

So far, we have considered stationary states for which the wall potential was left floating. 
We will now examine the effect of biasing the wall to a hxed potential below (ie, more 
negative than) the floating value (ffioat- Strictly speaking, the behaviour of the system in 
this case is not governed anymore by the ambipolarity condition at the wall, which was 
at the basis of the bounds obtained in Sect. |TT| However, the ambipolarity condition can 
be reintroduced using the fact that the ion current density is the same in both situations, 
because it is hxed by the boundary condition at the CS entrance. 
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(a) 


Figure 11: Collisionless case for Z1+, a = 2°, 
(j)^ below the floating value (fYioat ~ “2.5Teo. 
perpendicular to the wall. 



(b) 

and Uci = O.OSwpi with prescribed wall potential 
(a) Charge density profiles; (b) Ion flow velocity 


Still considering the electrons as perfectly magnetized np to the wall, we have 

{Vx) float'^i,float 


(Vx) 


bias 


'^i,bias 

= (n||)esina; 


n 


w 

e,bias 


n 


e,float 


( 11 ) 


‘^e^bias pbiasj ^ \ ^e,bit 


bias 


leading to the modihed bonnd 


1 • 1 / \Wi 

^(4^ float 4^bias^ 

= sm a (n||)g exp 

o 

_1 


( 12 ) 


Unsnrprisingly, the bonnd on the ion flow velocity becomes less and less restrictive as the wall 
potential is set to lower values. For a given target velocity, the corresponding critical angle 
decreases accordingly. Starting from a floating case, with a given (small) angle a for which 
the DS has nearly vanished, we can expect it to reappear as ((thias is decreased. Considering 
for instance the deuterium case with a = 2°, for which ~ —2.5Teo, several biased-wall 

simulations were performed with different values of 0^^. An increase of the charge density 
near the wall is indeed observed (Fig. [TT^), resulting in the growth of the electric held (not 
shown here) and the ion how velocity perpendicular to the wall (Fig. [IHd). 


We also analysed the case of a wall biased at a potential above (ie, less negative than) 
the boating value, a situation relevant to tokamak divertors where the divertor tiles may 
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(a) 


(b) 


Figure 12: Collisionless case for D^, a = 2°, 
(j)^ above the floating value 4>Yioat ~ ~2.5Teo. 
perpendicular to the wall. 


and LOci = O.OSwpi with prescribed wall potential 
(a) Charge density profiles; (b) Ion flow velocity 


be biased positively with respect to the floating potential. The results are depicted in Fig. 


12 For grazing incidence {a = 2°), a small bias above the floating potential is sufficient 


to remove completely the charge separation near the wall or even to reverse its sign. At 
the same time, the ion velocity at the DS entrance drops well below the sound speed. The 
conclusion here is that, for grazing incidence, a small bias above the floating potential can 
remove any remnants of the DS, so that the transition to the wall is truly charge-neutral and 
subsonic. For almost normal incidence, the necessary bias would have been much larger. 


V. COMPARISON BETWEEN THE FLUID MODEL AND KINETIC SIMULA¬ 
TIONS 

The results of Stangeby in were obtained using simple fluid model that had been pro¬ 
posed earlier by Chodura [2] and Riemann [7j, and further developed in [8]. Although its 
predictions are basically correct, most notably the disappearance of the DS for low incidence 
angles, it would be interesting to test its limitations by comparing the fluid results to those 
of our kinetic code. 

Taking the velocity moments of Eq. (§ up to hrst order yields the following fluid system 
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in the stationary state 


dxiuiU^) = - Uio) 

dx Pxx 


UxdrUr = 




UxdxUz = 


C/3^0 UJytlz U2 

mirii rrii rii 

dxPxy 1 ^iO 
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miTii j Ui 

9xPxz 1 , ^iO 


mifii 


OJyVjx ^x^y 


-Uy 


Ui 


(13a) 

(13b) 

(13c) 

(13d) 


where Uk = {vk),k = x,y,z, Ux = cjcisina, Uy = Udcosa, and rijo is the bulk density. In 
the Chodura-Rieniann-Stangeby (CRS) model for the collisionless magnetic presheath, we 
have r'i = 0 and two assumptions are made: (i) the non-diagonal components of the kinetic 


pressure tensor [terms in braces in Eqs. (13)c-d] are neglected and (ii) the xx component of 
the pressure tensor is assumed to follow an isothermal closure Pxx = To^i, with constant Tq. 
Combined with the quasi-neutrality relation and the Boltzmann law for the electron density, 


the system of Eqs. (13)a-d can be integrated easily |2l [10]. In Ref. [lO], the system is 
integrated in x starting from the CS exit. In our case, as the kinetic simulation encompasses 
both the CS and the DS, dehning the CS exit would require setting a somewhat arbitrary 
threshold on the charge separation. Thus, in order to compare our simulation results with 
the CRS fluid model, we integrate the fluid equations starting from the CS entrance at 
X = L. In order to compare with the kinetic results, the constant temperature Tq of the 
fluid model is set equal to the value of Txx = Pxx/n obtained from the ion distribution /*” 
at the CS entrance, given in Eq. (i‘. For clarity, as our notation differs from that used in 
Ref. [To], the explicit form of the CRS fluid equations is given in Appendix [Aj 


In Fig. 13 we compare the average velocity {vx) extracted from the simulation data with 
Ux computed from the fluid model for a few values of a. While the agreement is quite 
good for 0 = 3° and 5°, discrepancies appear for larger angles. It is important to note that 
those discrepancies arise before charge separation becomes signihcant, ie, when the plasma 
can still be reasonably considered as quasi-neutral {x > 10A/)). Proceeding to the same 


comparison for the y and z components of the average velocity (Figs. 14r and [l^), we 
observe quite similar discrepancies on but far larger and systematic ones for Uy on nearly 
the whole domain. Thus, as far as only the Ux prohles are concerned (and consequently also 


In our case, we have Tq = Tio(l — 0.55 sin^ a), where Tio = T^o is the parameter appearing in Eq. loL 
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2.5 



Figure 13: Ion mean flow perpendicular to the wall from the collisionless kinetic simulations (con¬ 
tinuous lines) and the CRS model (dashed lines), for various values of a, and deuterium ions. Note 
that the CRS model, being quasi-neutral, is not meaningful in the DS. 


the potential profiles), the predictions of the fluid model in the CS can be considered as 
rather good for the lowest range of incidence angles. The somewhat large and systematic 
discrepancies observed for the other velocity components would require closer scrutiny. They 
probably arise from the violation of both assumptions made in the fluid model. 


To rehne our comparison, we computed, from the kinetic simulations, the various terms 


entering the y and 2 ; components of the momentum balance equations (13)c-d. The com¬ 
parison indicates that the contribution of the non-diagonal terms of the pressure tensor is 


not negligible. Focusing in particular on the equation for Uy, Fig. 15 shows that the term 
containing P^y is comparable to the other terms, even in the CS. In contrast, the term 
(not shown here) is indeed negligible. We emphasize the fact that the non-diagonal nature 
of the pressure tensor is not an artifact due to the choice of coordinates, which could be 
eliminated by using a held-aligned basis: although the distribution at the CS entrance is 
indeed separable in (uy, v^), this separability is lost during the transition. 


Let us now consider the validity of the isothermal closure for the Pxx component of the 
pressure tensor. In the normal incidence case, for which only the DS exists, the temperature 
Txx (ie, the variance along Vx) decreases as the ion population is accelerated towards the 
wall by the electric held. This well-known ’’acceleration cooling” [THl [IS] persists in the 
magnetized case. More importantly, as the electric held prohle is spread out with decreasing 
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(a) 


(b) 


Figure 14: Same as Fig. 13 for Uy (a) and Uz (b). Note that, in the left panel (a), the continnons 
curves for a = 3° and 5° are virtually indistinguishable. 



Figure 15: Various terms of the momentum balance eqnation along y [Eq. 
the collisionless kinetic simulations, for a = 3° and deuterium ions. 


(13)c] compnted from 


a, Txx has a non-negligible variation over both the DS and CS. This is clearly visible in Fig. 
[I6| showing the evolntion of T^x relative to its valne at the CS entrance (we recall here that 
T^x^^ depends on a, see Sect. IVA). As a consequence, though the isothermal closure may 
be considered a reasonable approximation (outside the DS) for the large-to-intermediate 
angle range, it becomes clearly invalid in a large part of the transition layer for smaller 
angles of incidence. 

Having established that the isothermal closure does not £t the actual behavior of the 
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distribution for low a, one may hope to fit a slightly more general polytropic closure 
d{lnPxx) = 'jd{lnn). A typical constant polytropic coefficient 7 can be obtained by lin¬ 
ear regression for each value of a (Fig. mi- We observe a large variation with a, as can 
be expected when going from the two-scale behaviour at large a to the smoother transition 


at low a (see Fig. 16). For a = 90°, the CS disappears altogether and the quasi-neutral 
fluid model cannot be meaningfully compared to the kinetic results. Alternatively, one 
could compute a local polytropic coefficient 7 ( 0 ;) = d(lnP)/(i(lnn) [ 20 ], but this yields very 
large variations over the domain and with a, and is prone to numerical instability in the 
low-gradient zones. 


We also tried to use the computed exponent 7 to improve the match between the kinetic 
and the fluid models (using, in the latter, a polytropic equation of state, P^x 'n[), but this 


does not seem to work well for Ux (Fig. 18). The prohles of the mean velocities along y and 
z are not improved either, which is not surprising as their discrepancy with the kinetic data 
comes primarily from the assumption of isotropic pressure. The main conclusion here is that 
it is not possible to match the kinetic simulation data with a simple polytropic closure. 


All in all, the comparison between the simulation results and the predictions of the fluid 
model leads us to conclude that: (i) a rather good agreement is obtained for the Ux prohle 
(and consequently for the potential prohle) for the lowest values of a, but (ii) a worse 
agreement is observed for the other components of the mean velocity, due to the violation 
of some of the assumptions of the huid model. 


VI. COLLISIONAL SIMULATIONS 

In the preceding collisionless simulations, the held-aligned ion how velocity at the CS 
entrance was imposed through an ad-hoc boundary condition. In order to ensure that 
such results are not specihc to the collisionless regime, we performed a series of collisional 
simulations. In this case, the simulation domain is much larger (typically 2 x 10 ^Ad) in 
order to encompass the full transition from the the plasma bulk to the wall. The ion velocity 
distribution in the bulk is an isotropic Maxwellian with temperature Tjo = Tgo. Then, the 
distribution function at the CS entrance is no longer imposed as a boundary condition, 
but rather arises self-consistently in the collisional presheath located upstream the DS. A 
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Figure 16: Spatial variation of the temperature T^x normalized to its value at the CS entrance 
Txx^j for collisionless simulations with deuterium ions. 




(a) (b) 

Figure 17: (a) Determination of an average polytropic index 7 by linear regression of InP^^^; = 
/(ln(n)) for a G {2°, 30°, 60°}. Simulation data are plotted as dashed lines and regression results 
as continuous lines; (b) Behaviour of the average index 7 with a. 


thorough characterization of the transition, using the same kinetic model, was performed by 
Devaux et ah [21]. Here, we will focus on the question whether collisions modify the results 
obtained in the collisionless regime for grazing incidence. 

As in the collisionless simulations, parametric scans in a were performed for the same 
range of angles for ojd = 0.05upi, with three valnes of the collision freqnency Ui G 
{5 X 10“^, 10“^, 5 X 10“^} For this range of parameters the transition is charac- 
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Figure 18: Ion mean velocity profile {vx) for the kinetic simulation, the fluid CRS isothermal model 
(7 = 1), and the fluid polytropic model (7 = 1.74). Deuterium ions with incidence angle a = 3°. 


terized by the scaling Xd ^ Pi ^ Xm/p, where A^/p = Vth/^i- This is the intermediate Bq 
regime described in Refs. [m ED 122], for which the collisional presheath, Chodura sheath, 
and Debye sheath are well separated. 


As a preliminary benchmark, we use the collisional simulations to check the validity of the 
boundary condition that we prescribed at the entrance of the CS in the collisionless runs [Eq. 
(|^]. For this, we need a criterion to define the CS entrance. In the collisional presheath, 
the mean ion velocity is aligned with the magnetic held: the CS entrance corresponds to the 
location where this alignment breaks down. As a quantitative criterion, we took a deviation 
of 0.25° with respect to the angle of incidence a. The computed distribution functions are 


shown in Fig. [^and look very much like the prescribed distributions used in the collisionless 
runs (Fig. |^). 

We can now verify the robustness of Stangeby’s result in the collisional regime. First and 
foremost, we still observe a decrease of the charge density near the wall for decreasing angles 


of incidence (Fig. 20), with similar consequences on the electric held and potential prohles 
near the wall (not shown). The principal ehect analysed in this work is thus not destroyed 
by the presence of collisions. 

Second, the nearly linear dependency of the wall charge density with sin a (which was 
observed in the collisionless case, see Fig. 1^) is slightly perturbed by the collision terms as 


shown in Fig. 21 (note that here the charge density is normalized to the value hq in the bulk 
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^x^thi 


Figure 19: Ion velocity distribution functions at the CS entrance, for various angles of incidence, 
t'i = 5 X and deuterium ions. The position of the CS entrance is indicated in the inset. 

In order to facilitate the comparison, each distribution function is normalized in such a way that 
it has the same maximum as the prescribed collisionless distribution of same incidence angle (Fig. 

[^). 


plasma, whereas in the preceding sections the normalization value Uq referred to the density 
at the CS entrance). A marginal sign inversion of p near the wall can even be observed in 
the (a = 2°,z/j = 5 X 10~^Vth^D^) case. Despite this perturbation, the ion perpendicular 
flow as a function of a may still be roughly htted by the same semi-empirical law as in the 


collisionless case (Fig. 22). 

Last, let us extend the analysis of the various terms entering the fluid momentum balance 


in Eqs. (13)b-d. Setting aside the additional impact of the friction terms specihc to our 
collision model, we still observe a non-negligible impact of the non-diagonal term of the 


pressure tensor P^y in the fluid momentum balance along the y axis (Fig. 23). As was the 
case for the collisionless regime, the P^z cross-term (not shown here) is indeed small outside 
the space-charge region near the wall. 
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(a) (b) 

Figure 20: Charge density profiles near the wall for two collisional simulations with deuterium ions. 
The collision frequencies are i/i = 5 x 10~‘^vth/ (a) and z/j = 5 x vth/ (b). 



Figure 21: Charge density on the wall as a function of the incidence angle a, for three values of 
the collision frequency r'j. Deuterium ions. 


VII. CONCLUSIONS AND PENDING ISSUES 

The main focus of this paper was on the observation, made by Stangeby ng, that the Debye 
sheath should disappear when the plasma is immersed in a magnetic held with grazing angle 
of incidence with respect to the wall. Stangeby’s result was deduced from a simple ID huid 
model with Boltzmann electrons and isothermal closure for the ions. Thus, it was worth 
to check whether the result holds under less stringent conditions on the ion model, namely 
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Figure 22: Ion flow velocity on the wall as a function of the incidence angle a, for a collisional case 
with frequency i/j = 5 x lQ~^vth/^ d- The numerical coefficient k = 7.31 of the semi-empirical law 


Eq. (10) is computed by exact interpolation from the case a = 90°. Deuterium ions. 



Figure 23: Various terms in the momentum balance equation along y, Eq. 
simulation with D'^ ions, a = 3°, and i/j = 5 x ^d- 


(13)c. 


Collisional 


using a kinetic rather than fluid approach. 

Our calculations showed clearly that the main result holds: the charge separation pro¬ 
gressively disappears for smaller and smaller angles of incidence, and the ion flow velocity 
perpendicular to the wall is limited to subsonic speeds. Though no critical angle arises due 
to the lack of singularity at the DS entrance in the kinetic model, the overall behaviour is 
consistent with the predictions of Ref. [lO]. We also conhrmed the increased spreading, with 
decreasing a, of the electric held and plasma density over distances of several Larmor radii 
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from the wall. These features appear in both collisionless and collisional simulations, and 
may thus be considered as robust, provided the scaling Xcou 3> p* 3> A^) is satished. 

As noted by Stangeby [10], the spreading of the electric held and plasma density further 
from the wall (compared to what is usually expected from simpler models) has important 
consequences on the recycling of sputtered particles in a tokamak edge. It should be taken 
into account, whenever possible, in the computational codes that deal with plasma edge 
recycling. 

Further, by comparing the kinetic and huid prohles, we found that, although a rather good 
quantitative agreement on the ion how velocity perpendicular to the wall can be obtained 
for small angles, the assumptions of a scalar pressure tensor and isothermal closure in the 
huid model are clearly violated. These hndings point at the limitation of the huid models 
usually employed to study this type of scenarios. 

Finally, in all simulations apart from the most collisional ones, we observed a rather robust 
linear scaling of the charge density at the wall with sin a. As a consequence, the value of 
the ion mean how velocity perpendicular to the wall obeys the simple semi-empirical law: 
{vx)Y = 'i^t/ie/\/^sina;/(l -|- Ksina), where k is a coefficient that can be determined from 
a single simulation at normal incidence. 

All the previous considerations are correct as far as the various simplifying assumptions 
made both in the huid and kinetic models are satished. The hrst concerns the electrons, 
which were assumed to be perfectly magnetized up to the wall and to follow a Boltzmann 
law. For very small angles of incidence {a < 1°), these assumptions cease to be valid and 
the electron dynamics should be treated with a fully kinetic model. 

A second assumption lies in the reduction of the system to one dimension in space. For 
divertor targets, the determination of the CS and DS structure near the inter-tile gaps would 
require at the very least a two-dimensional model in space, encompassing the full incidence 
plane of the magnetic held [ie, the plane {x,y) in our geometry, see Fig. in order to 
properly determine both the structure of the electric held and the particle hows in those 
regions. Of course, an extra spatial dimension would increase dramatically the complexity of 
the present kinetic code. Nevertheless, it is an important feature that needs to be addressed 
for quantitative comparisons with tokamak measurements. 
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Appendix A: Collisionless fluid model 


The following relations are established from the fluid system (13) in the collisionless case 


iyi = 0) using a diagonal pressure tensor {P^y = Pxz = 0) and an isothermal closure P = nTo. 
Exact neutrality rii = rig and a Boltzmann law for electrons are assumed. The integration 
of the system follows the same pattern as in Refs. mm, the only difference being in the 
fact that no assumptions were made on the value of the boundary conditions (ie, they are a 
priori unrelated to Cg). 

Starting from a reference point Xq with fluid velocities {u^o < 0,Uyo,Uzo), the position 
Xi < Xq where reaches the value u^i is obtained through the integral expression: 


'^xl 


Uci cosQ:(a;i — Xq) = — 


u{l- du 


UxO 
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and Uq = uIq + u^q + u^q, cl = (Tq + Tgoj/mj. The above relations are obviously valid only 
as long as D{u) does not vanish in the integration range. Here, the Cg factor arises solely 
from the isothermal closure for the ions, and does not depend on the boundary conditions. 

From a numerical point of view, the velocity prohle Ux{x) is reconstructed as follows: 
a uniform discrete velocity grid (m„, n = 0... A) is generated between uq = and 

mtv = max(— —Ubound), where Uging is the singular velocity cancelling D{u) and Ubound 
is the velocity bound obtained from Eq. Q. Starting from [uo,ui] Eq. (Al) is integrated 
over each pair [un, Wn+i]- The end result is a sequence [xq, ..., xn] of positions matching the 







velocities [mq, ... ,un]- The Uy profile is obtained directly using Eq. (A3). The velocity Uz 
is recovered from using 


and the electrostatic potential 
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